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Abs TRAC I 


To control a lightweight manipulator, a validated dynamic model based on the 
Equivalent Rigid Link System (ERLS), is used to formulate a control algorithm. The 
required torque to drive the manipulator to a desired angle is calculated. Since the 
control system includes an electrohydraulic actuator, the required current is also 
calculated by the inverse dynamics of the hydraulic servo. A computer simulation is 


performed with various loading conditions, gain values, and desired angles. 
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I. INTRODUCTION 


A. MOTIVATION 


The robotic manipulator plays an important role in a wide range of applications. 
Present and possible future applications focus on implementing automation in 
environments which are hostile to humankind. More demand 1s being placed on robots 
with higher speeds, smaller actuators, and higher payload capacity in applications. 

In the past, efforts were made to control manipulators using rigid body 
assumptions. However, the demands for lightweight and higher performance 
manipulators have led to the consideration of flexible- body models. A new controller 


is needed for the flexible manipulators. 


B. NATURE OF THE PROBLEM 

Control of most manipulators is based upon dynamic models assuming perfectly 
rigid links. Accordingly, accurate manipulator motion control is dependent upon the 
great stiffness of the structural members. This assumption does not always hold, 
particularly at high-speed and heavy-load operations where a linkage may undergo 
severe elastic deformation due to the inertia. It can no longer perform its intended 
function in actual operations. 

Lightweight flexible manipulators promise the advantages of increased speed of 
Operation, ease of transportability, reduced material requirement, reduced power 
consumption, lowered mounting strength and mgidity requirement, and lower overall 
cost [Ref. 1]. On the other hand, the oscillatory behavior of the flexible links makes 
the end-point motions difficult to be controled. The central problem lies in the control 
of the end-point motion by applying the proper torques at the actuating point. 
Flexibility effect on manipulators is a necessary consideration in the move toward 
designing lighter, faster and more accurate robot systems. A dynamics model of 
manipulators with flexibility is needed for design and control purposes. To benefit 
from the advantages of lightweight, flexible manipulators, it is also necessary to 
implement a control design capable of achieving end effector positioning accuracy and 


Stable control. 


C. LITERATURE REVIEW >. 

The earliest studies concerning the control of flexible manipulators began in the 
early to mid 1970's. There has been considerable research recently in the development 
of flexible manipulator control strategies. A brief survey of this research follows. 

Book [Ref. 2] introduced a modeling approach utilizing 4 x 4 transformation for 
control purposes. His model was limited in the assumption that the mass of the 
manipulator is negligible compared to the mass of the load. The assumption is 
acceptable in space applications where the manipulator is very light and operates at 
low speeds, but the assumption can not be applied to the industrial robotic 


manipulators in most cases. 


Sunada and Dubowsky [Ref. 3] applied 4 = 4 transformation matrix and the 
Lagrangian approach to model the flexible system with geometrically irregular links by 
finite element techniques. Component Mode Synthesis was used to reduce the number 
of generalized coordinates and improve the computational efficiency without losing 
much accuracy. Based on this scheme, they can obtain the flexibility and mass 
properties of the links through available finite element programs. But, this model 
ignored the effect of the small motion interaction on the large motion. The 
consideration of the coupling effect between the small motion and the large motion 
was not completed. 

Chang [Ref. 4] introduced the Equivalent Rigid Link System(ERLS) dynamic 
model of flexible manipulator. Global motions were separated into large motions and 
small motions. The ERLS described the kinematics of the large motion. The small 
motion displacements were described relative to the ERLS. The finite element method 
was used to discretize deformations and Lagrangian formation was used to derive the 
equations of motion. Two sets of coupled, non-linear ordinary differential equations - 
large motion equations and small motion equations resulted. The set of large motion 
equations were non-linear in both the large and small motion variables. The set of 
small motion equations were linear in the small motion variable but non-linear in the 
large motion variable. The model is complete being able to describe large motions, 
small motions and their coupling. The model presented the computational potential by 
using the Sequential Integration Method developed. 

Petroka [Ref. 5] performed an experimental validation of the ERLS model. He 
designed a hydraulically actuated, single - link flexible arm. A cubic shape function 
was assumed to present the transverse displacement of the flexible manipulator. Tip 


position was determined from the motion picture studies. 


1] 


Gannon [Ref. 6] served a continuation of the experimental validation of the 
ERLS model. A piezo-resistive accelerometer and strain gages were used for tip 
position and strain information. The natural mode shape functions of a beam was used 
to present the flexural motion of the single-link flexible arm in the Gannon’s study. 
Results of the validation suggested the potential usefulness of the model in determining 
t1p position. 

Canon and Schmitz (Ref. 7] reported on experimental studies of the “End - Point 
Control of a Flexible One Link Robot”. Their controller design was based upon a 
linearized and truncated model. Measurements available for control included rigid body 
angle and velocity, strain measurements and direct optical measurement of end-point 
position. Controller design was discussed from a variety of standpoints using both 
classical and modern control methods. The authors reported that the non 
-colocated(end -point position) feedback improves the speed of response at the cost of 
reducing the stability margins. 

Flexibility effect is a very important consideration in designing lighter, faster and 
more accurate robot systems. In this research, the ERLS dynamic model is used to 
design a control svstem since the model presents a complete consideration of the 
system motion and exhibits computational superionty. Nonlinear motion is 


compensated in the control algorithm. 


D. RESEARCH OBJECTIVE 

The objective of this research is to design a controller for single-link flexible 
manipulator. The validated dynamic model-the Equivalent Rigid Link System - is used 
to formulate the control algorithm. The control algorithm involves the computation 
of the required torque to drive the manipulator achieving a desired angle of the end of 
the flexible manipulator. Since the control system includes an electrohydraulic actuator, 
the required current is also calculated by inverse hydraulic dynamics. The ways to 


reduce the actuator effort and structural vibrations will be discussed. 
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Il. MODEL FORMULATION FOR A SINGLE LINK FLEXIBLE 
MANIPULATOR 


A. PHYSICAL PLANT 

A flexible arm was designed by Petroka (Ref. 5] for the validation of the 
Equivalent Rigid Link System Model. The flexible arm shown in Figures 2.1 - 2.2 is a 
one meter long flexible structure which can be bend freely in the vertical plane but is 
stiff in the torsion and horizontal bending. The basic configuration of the arm includes 
two parallel, flexible steel, flat bars connected by the thin steel strips to transverse steel 
bridges. External loading is attached to the arm tip end transverse bridge by the 
securing of the load on the four welded studs shown in Figure 2.3. Table 1 shows the 
geometric and mass properties of the flexible arm. The flexible arm is driven by a 
electohydraulic system shown in Figure 2.1. The system consists of a York hydraulic 
power unit, a Berd- Johnson 3 - axis Hyd-Ro- Wrist, a Moog 760-100 servovalve, and 


associated Moog servocontroller and high pressure filter assembly [Ref. 5]. 


B. EQUIVALENT RIGID LINK SYSTEM FOR A _ SINGLE-LINK 
MANIPULATOR 


The basic idea of ERLS [Ref. 4] 1s to separate the motion of the flexible link 
system into the large rigid-body motion and the small-motion displacement due to the 
flexibility of the structure. The large motion is defined by the movement of ERLS for 
a single - link flexible manipulator. The small-motion displacements are described with 
respect the ERLS. A schematic diagram of the ERLS for a single - link manipulator 1s 
shown in Figure 2.4. 

The three generalized coordinates chosen are the large motion joint variable,8, 
and two nodal displacements,V(0)and (0) ,which can be measured at the end of the 
link. The large motion joint variable ,8, is the angle between the ERLS and the 
horizontal axis of inertia coordinate. Homogenous transformations is applied for a 
deformed point on the single-link to obtain the absolute position of the point. The 
displacements for each point along the flexible arm are the functions of the location 
and time; therefore, it is necessary to discretize the deformations in terms of generalized 
coordinates in oder to apply Lagrangian’s equation. The techniques of the Finite 
Element Method (FEM) are utilized for this discretization. The nodal displacements at 
the end of the single-link represent the small motion. 
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ARM PARAMETERS 


| Arm Length (L) 0.9985 


Beam Thickness(t) 0.003175 m 
PII peuicen 0.05715 m 
Modulus of Elasticity(E) pe) N/m? 


Density/unit Volume( p) 7861.05 kg'm? 


After describing the kinematics relationships between the large and small 
motions,kinetics is introduced to complete the derivation of the equation motion. 

The Lagrangian dynamics approach is used for the derivation due to its 
straightforward and systematic nature which is helpful in the analysis of complex 
systems. The kinetic energy is the sum of the kinetic energies of each link,actuator,and 
any loading. The potential energy is contributed the elastic strain energy and from 
gravity. Generalized forces are included due to any applied forces and damping forces. 
After making considerable effort in mathematical manipulations, rearrangements, and 
simplifications, the Lagrangian equations vield two sets of non-linear, coupled, second 
-order, ordinary differential equations. One set describes the large motions and the 
other set describes the small motions. 

These equations are represented as follow; 

Mgq® + MgnU =F 


+TORQUE (2.1) 


q q 
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Mng@ +MpnU + K,U =F, . (2.2) 


where 

Mag = 1 X 1 effective inertia matrix for large motions 

Mon = | X 2 coupled inertia matrix of the small motion effect 
on large motions 

Ming = 2 X 1 coupled inertia matrix of the large motion effect 
on the small motions 

Mnn = 2 * 2 effective inertia matrix for small motions 

Kn = 2 X 2 stiffness matrix for small motions 

Fg = | X | load vector for large motions 

FL, = 2 X 1 load vector for small motions 

0 = generalized coordinate of the large motions 

U = 2 X 1 generalized coordinate of the small motions 


The vector of small motion is limited to the transverse displacement of the end of 
the link. Axial deformation and torsion are neglected in this research for single-link 


manipulators. We assume that the displacement of the end of the link is small. 
Tan!(V/L) ~ V/L (2.3) 


where V is a shorthand of V(0). 

We also assume that the slope of the arm tip is not measured for the simple 
feedback control system. The tip displacement V(0) is the only generalized coordinate 
chosen for small motions. Then all coefficient of the matrixs of non-linear, coupled, 
second-order ordinary differential equations is reduced to | X 1 matrix. A total angle 
of the end of the link can be represented as follows; 


@M=6+ V/L (2.4) 


® = 0+ VIL (2.5) 


ig 


® = § + VIL s (2.6) 


A more detailed derivation of the equations of motion 1s included in Appendix A. 


C. SHAPE FUNCTION DERIVATION 

The natural- mode shape function of a beam is needed to present the flexural 
motion of the single-link flexible arm [Ref. 4]. The flexible manipulator arm is modeled 
as a continuous Euler-Bernoulli cantilever beam. The transverse displacement for the 


single link flexible arm is represented in the following forms; 


Ulx,t} = ay (t) xy (xd tay (1) XW) (2.7) 
= ay, {A,’ {cosB)x + cosh Bx} + {sinB,x + sinhB,x}} 
+ ay [A ) {cosB4x+cosh B5x} + {sin Bx + sinh B5x}} 


where 
BL = 1.875104069 
BL = 4.694091133 
sin B-L + sinh B-L 
= S8BiL + sinh BL aa 
Cos) Pl iepcoshah.k 


Upon reduction, Gannon [Ref. 6] yielded a 3 X 2 shape function matrix ¥, and 
the 2 X 1 nodal displacement vector, U. a 3 X 1 deformation vector, d , was 


expressed as follow; 


d=WU = 0 =( 0 0 V(0) 
0 0 0 (0) 
V(x) M N 


The coefficents Mand N are defined as; 


M = {Cy {Ay’ ( cosB,)x+coshB,x }+{ sinB)x+ sinhB ,x}} 
+ Cy { Ao’{cosB4x+ coshB5x}+{sinB, x+sinhB>x}}} 
N = es {Ay { cosB,)x+ coshB ;x ae sinB)x+ sinhB ,x}} 
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+ Cy { Ay'{cosBax+ coshB>x}+{sinB, x+sinhB,x}}} 
where 


ee 
| 


= 2 B,/{4 Ay’ B,-4 A,’ By 
C> =-2 8 ,/{4 Jay p> -4 Ay By, 
C3 =-2 Ayirq Ay Bg - 4 Ad’ By 
Cy = 2 ALif{4 oa B5 - 4 NDS By 


In this research, since the slope function of the arm tip is not measured, the 
shape function, VY , is reduced toa 3 X 1 matrix. Now the 3 X | deformation vector 
can be rewritten as follows; 


d=¥U= 7 0 )=f 0 [% ] 
0 0 
V(x) M 


New shape function matrix is now in a convenient form for computer coding. A 


more detailed derivation of the shape function was included in Gannon’s work. 


D. HYDRAULIC ACTUATION 

The single-link flexible manipulator uses electrohydraulic actuation to drive the 
plant. The inclusion of the electrohvdraulic power system involves the transformation 
of an input current to an output torque. Figure 2.5 illustrates the relationship of 
hydraulic dynamics to the overall system. 

The hydraulic dynamics is represented by servovalve dynamics and actuator 
dvnamics. Moog, the manufacture of the servovalve, simplified the description of 


servovalve dynamics to a single equation [Ref. 8]. 


O=1K/ Pe, (2.9) 
where 
Q = flow delivered from servovalve 
[ = input current 
K = valve sizing constant,contributes to hydraulic system damping 
P., = valve pressure drop,P.-Py 
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Actuator dvnamics consists of a form of the continuity equation and the torque 
output equation. Actuator dvnamics consists of a form of the continuity equation and 


the torque output equation [Ref. 9}. 


0 = D260] Guana suse (2.10) 
m ine we a 6. 
Ta = Tet ee : (2.11) 
W here 
Q = flow delivered to actuator 
De 8 = flow causing actuator rotation 
B. = effective bulk modulus of fluid 
. = total compressed volume including actuator lines and chambers 
ve Py i (4 Be ) = compressability flow 
Tg = torque required to overcome inertfa and move the load 
nN, = torque efficiency 
Fe = load pressure drop 
Da, = motor displacement 
CemPyp = leaking flow in actuator 


tJ 
tJ 


The preceding hydraulic dynamic equations and relationships are incorporated 


into the main computer program for solving the dynamic equation of motion. 


ZS 


lil. CONTROL SCHEME 


A. REQUIRED TORQUE 

In feedback control system, desired angles are compared with the measurements 
of the actual system responses. The differance is used via the chosen control 
algorithm. The control algorithm calculates the required torque to derive the 
manipulator achieving a desired angle of the end of the arm. 

Recall the equations that can be written as two sets of equations - large motions 
and small motions equations - and the total angle equation of the end of the arm for 


derivation of the required torque equation. 


Mgq@ + MgnV =F + TORQUE (3.1) 
Mag ® +Man ¥ + Ky V = Fy (3.2) 
where 
M qq 7 1 X | effective inertia matrix for large motions 
M qn = 1 X 1 coupled inertia matrix of the small motion effect 
on large motions 
Mang = | X 1 coupled inertia matrix of the large motion effect 
on the small motions 
Man = 1 * I effective inertia matrix for small motions 
a = | X | stiffness matrix for small motions 
Fg = |] X | load vector for large motions 
F, = I *X 1 load vector for small motions 
9g = generalized coordinate of the large motions 
V = | X 1 generalized coordinate of the displacments 
86 =@M-V/L (3.3) 
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® - V/L 


- 
ll 


@=@-VIL 


(3.4) 


(3.5) 


Equations 3.3-3.5 can be substituted into Equations 3.1 and 3.2 and rearranged as 


follows; 


Mag (D-V/L)+ Man V-F. = TORQUE 


1 


oul) M,, Vt K, V-F, = 0 


After eliminating V, a governing equation for ® is given as; 


(Mgq -D Mag) ®+(D Fy -D Ky- Fy) V = TORQUE 


where 


5 o Man Magg /L) 


vie (Mig /L)) 


Let 


N = Mog -D Mac 


Equation 3.8 can be expressed in short as follows; 
N@® + FC(V,V,0,6) = TORQUE 


where 


BORG, 0) = fey 0) Ko ey 
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(3.6) 


(3.7) 


(3.8) 


(3.9) 


(3.10) 


@ is measured and fed back to controller.. The required torque to drive the 
manipulator achieving a desired angle of the end of the arm is calculated via a 
controller. [he inverse problem is applied to calculate the required torque. Figure 3-1 
shows a schematic diagram of the feedback control system. 


The equation of the required torque is represented as follows; 


ds = N( @®gt i (®, -D) + Ky (Dz -@))+ EC Gas 
where 
T, = a required torque 


®, = a desired angular acceleration 
® = a total angular acceleration 
db, = a desired angular velocity 

® = a total angular velocity 

@, = a desired angle 

@M= a total angle 


Ky = a velocity seam 
K 


p74 position gain 


The torque of the Equation 3.10 is equals to the required torque of the Equation 


3.11. System errors of a second order control system can be represented as follows; 


(O,-®) + K,(®, -®) + K, (®g -®) = 0 (3.12) 


B. REQUIRED CURRENT 

Since the control system includes an electrohydraulic actuator, the inverse 
hydraulic dynamic involves the transformation of input torque to an output current. 
For the input current to the hydraulic dynamics, the inverse hydraulic dynamics can be 
applied to calculate the required current according to the required torque. A desired 


load pressure drop and desired flow should be calculated by the inverse hydraulic 
dynamic equations. 
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A Schematic Diagram of the Feedback Control Svstem. 


Figure 3.1 


Q, = D8 + Cun Pre t(V; Pr, 4 Be) -- (3.14) 


1 = Qqg (KVP, ) (3.15) 
where 
T, = a required torque 
P;, = a required load pressure drop 
Q, = arequired flow delivered to actuator 
I. = a required current 
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IV. RESULT 


A control algorithm for a single-link flexible manipulator is based upon the 
Equivalent Rigid Link System Model. The control algorithm calculates the required 
torque to control the end position of the manipulator. Comparisons of the required 
torque and required current are made in terms of the three parameters (i.e., loading 
condition, gain value, and desired angle). Additionally, the way to reduce the required 
Currents and required torques is discussed. The effect of the current saturation to 


system dynamics behavior is also presented. 


A. INPUTS AND MOTION OUTPUTS 

The loading conditions are no load, a 2.115 kg ,and a 4.233kg load. The 
external loading 1s attached to the arm tp end transverse bridge by securing the load 
on the four welded studs. The initial condition for all runs is the horizontal position of 
the flexible arm. The initial tip deformations as determined by the Euler-Beroulli 
cantilever beam theory are converted to an angle. These angles in Table 2 are the 


initial conditions input into the simulation program. 


TABLE 2 
INITIAL CONDITIONS 


Load ( Kg ) Angle (radian) 


2.115 -0.14628 


4.233 -0.23155 


Only point to point control in open space is considered. An input variable to the 
closed-loop system is the desired angle of the end-point. The desired angular 


acceleration and the desired angular velocity are zero in this research. To compare the 
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system responses, two step inputs, Which are two desired angles of the end-point,(i.e., 1 
radian and 1.5 radian), are applied respectively. 

The outputs of the system responses for the step input of the desired angle are 
total-angle accelerations, total-angle velocities and total angles of the end-point, which 
are also fed back to the controller. First, the desired angle of 1 radian is input as the 
step function. Figures 4.1-4.3 show the movements of the end-point for all loading 
conditions. Each figure has three curves according to selected gain values for all 
loading conditions. In this research, a velocity gain ( K,) and a position gain (Ky) are 
selected for critically damped systems. In other words, the movements of the total 
angle of the end-point will not exceed the desired angle of the end-point. The position 
gain can be expressed in terms of the velocity gain for the critically damped systems as 
fellows; 


=K2 
Kem 


Table 3 shows the three cases of gain values for this research: 


TABLE 3 
GAIN VALUES 


. 
| 
[> 


For the evaluation of the system performances for the step input, a settling time 
corresponding to a 2 % tolerance is measured in terms of time constant( T = 1 / 
Keer Settling time( t, = 4 T ) are 4 sec, 2 sec and 1 sec according to the position 
gains 1,4 and 16. All results of the total angles of the end-point meet error criteria of 
2 % according to the gain values. The figures of the total angle of the end-point with 


larger gain values show faster movement. Another observation is that the total angle 
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of the end-point with the same gain value have the same movement in terms of the 
time constant for the different loading conditions. These results indicate that the 
nonlinear plant is controlled to perform a linear behavior. Errors of the system can be 


represented as the second-order control system as follows; 
(Dy-D) + Ky, (®y -&) + K, (®q -®) = 0 


Figures 4.4-4.5 show the total-angle velocity and total-angle acceleration 
approaching the desired angle of | radian for the 4.233 kg loading condition. All 
initial velocities start from zero and being increasing to achieve the desired angle with 
fast movement as gain values are increasing during transients and converge to zero 
velocity as the total angle of the end-point approachs to the desired angle. 

Figure 4.6 shows the comparison of the total angle response for the step inputs 
of | radian and 1.5 radian with the position gain of 4, the velocity gain of 4 and 2.115 
kg loading condition. The total angle movement of the 1.5 radian is faster than the 
movement of the desired angle of 1 radian. But both of the total angle response have 
same settling time since they have the same position (4) and velocity gain values(4). 
Figures 4.7-4.8 show the comparisons of the total-angle velocity and total-angle 
acceleration between the | radian and 1.5 radian of the desired angle with same 
position gain for the 2.115kg loading condition. The figures of the total-angle velocity 
and total-angle acceleration indicate that total-angle velocities and accelerations are 
increasing as desired angle is increasing with same gain value during the transients. 
The maximum of total-angle velocities and total-angle accelerations for the desired 


angles are shown in Table 4. 


B. REQUIRED TORQUES 

In the feedback control system, the desired angle of the end-point is compared 
with the measurements of the actual system responses. The differences are used via the 
chosen control algorithm to calculate the required torque. In other words, the control 
algorithm calculates the required torques to drive the manipulator to achieve the 
desired angle of the end-point according to the chosen gain values. The equation of 


the required torque was represented as follows; 


T. = N(@j + Ky (®q -®) + K, (gq -®))+ FC 


ol 


TABLE 4 


MAXIMUM OF TOTAL-ANGLE VELOCITY AND TOTAL-ANGLE 
ACCELERATION 


| Desired Maximum Angular Maximum an va 
| Angle(rad) tad sec) Acceleration(rad sec’ 
1.0 ar ee 974 Sel 


Figures 4.9-4.11 show the required torques to drive the manipulator in order to 
achieve the desired angle of | radian for all loading conditions. The desired angle of | 
radian 1s input as a step function. All figures have three curves according to the three 
gains in Table 1. Table 5 shows the maximum of the required torques according to the 
gain values and the loading condition. The figures for the required torques and Table 
5 show expected results; the required torques during the transient period are increasing 
to deliver the load and to overcome the inertia moment as gain values are increasing, 
and rapidly settle to the steady-state value. These required torques are related to the 
movements of total angle of the end-point. For the faster movements of the total 
angles, more required torque should be applied to the physical plant. Because of no 
steady-state error, the required torques at the steady-state ( the total angle reaches to 
the desired angle ) have same values for the same desired angle even though the gains 
are different. 

Figure 4.12 shows the required torques for desired angles(i.e.,1 radian and 1.5 
radian). The comparison for the desired angles indicates that the maximum of the 
required torque depend upon the movements of the total angle is increasing as the 
desired angle 1s increasing. 


C. CURRENT CALCULATIONS 
The control algorithm includes an electrohydraulic actuator. The input current 
to the manipulator can be calculated from the inverse electrohydraulic dynamics 


according to the required torque. Figures 4.13-4.15 represent the graphical results of 


Be 


ewe 


vee ei REQUIRED TORQUE 





a 


148.11 


] 


the required currents for all loading conditions. Each figure has the three curves 
according to three gain values. Figures indicate that the required currents have 
fluctuations based upon required torques during the transients, and settle to the 


Steady- state values. 


D. THE REDUCTIONS OF THE REQUIRED TORQUE 

The magnitudes of the required torques and required currents are related to the 
selection of size of the electrohydraulic actuator. The selections of the electrohydraulic 
actuator which depend upon the required torques and required currents, are also 
related to the weight,and cost. To reduce the required torques and required currents is 
very important in terms of the selection of the size of electrohydraulic actuator. In this 
research, a planned desired angle trajectory is used as the input replacing the step 


input. 
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The planned input is the combination of ramp and step function applied to the 
4.233kg load with the position gain of 16. The initial point of the ramp function 
coincides with the initial angle of the flexible arm. The slope of the ramp function is 
selected as 2. Figure 4.16 shows the planned input of the desired angle as the 
combination function. Figures 4.17-4.19 show comparisons of the results between the 
step input and planned input. The difference of the maximum required torque for 
input functions 1s tabulated in Table 6. 


| TABLE 6 
MAXIMUM TORQUE FOR INPUT FUNCTIONS 


Planned Input 106.58 


Table 6 indicates that the plan of the input function for the desired angle could 
results in the reduction of the maximum required torque and maximum required 
current. The total angle responses might be slower. We can expect different reductions 
of the required torques or required currents for same loading conditions and gain value 


when We change the input function of the desired angle. 


E. THE REDUCTION OF VIBRATION 
The system damping is considered to reduce the structural vibration. A damping 
ratio 1s applied to the small motion equation as follows; 


Mag ® + Man ¥ 26 nn oe a 


The no load fundamental frequency is predicted to be 2 HZ. A damping ratio( ¢) of 
0.5 is applied to no load condition with the position gain of 2. Figures 4.20-4.23 show 
graphical results of the comparison of the no damping case with the damping case that 
has a 0.5 damping ratio. These result show that the movement of the end-point of the 


damping case is faster than the movement of the no damping. The system damping 
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yields the noticeable reduction of the vibration of small motion. The required torques 


of the system damping case show less fluctuations. 


F. THE ELECTROHYDRAULIC SATURATION 

The previous results for the required torques and required currents are for an 
ideal condition. When the electrohydraulic actuator is saturated at 15 miulliamp 
current, all current exceeded 15 milliamp will be limited at 15 mulliamp current. 
Figures 4.24-4.26 show the effect of the saturation for the 4.233 kg loading condition 
with the position gain of 16. The total angle responses with the saturated current show 
slightly slower movements compared to the ideal movement during transients. The 
required torques are also reduced due to the saturated current during the initial 


transients. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

The purpose of this research is to design a controller for a single-link flexible 
manipulator using the Equivalent Rigid Link Svstem. The control algorithm 1s 
formulated calculating the required torque to control the desired angle of the end- 
position. Comparisons of the system responses and required torques were made in the 
different conditions. 

The results indicate that system responses and required torques are related to the 
selected gain values for all loading conditions. The movements of the end-point with 
higher gain values are faster to achieve the desired position. For different desired 
positions with the same gain, the end-point movements have the same pattern. These 
results indicate that non-linear motion is compensated in the control algorithm. 

The required torques are increasing to deliver the load and overcome the inertial 
forces during transient times as gain values are increasing. The system rapidly settles to 
the steady state. Because of no steady state error, the required torques at the steady 
state have same values for the same desired positions even though the gains are 
different. 

The maximum of the required torques and required currents are related to the 
capacity of the electrohydraulic actuator. The selection of size of electrohydraulic 
actuator affects the weight and cost of the control system. Since more demand 1s being 
placed on robot control with higher speeds and smaller actuators, it 1s important tc 
reduce the maximum required torques and required currents. The application of the 
planned input of the desired angle as a combination of ramp and step functions lead a 
noticeable reduction of the maximum required torques and required currents. However, 
slower movements of the end-point are vielded during the transient period. The 
maximum required torques and required currents can be adjusted by the slope of the 
ramp function. 

The results also indicate that the system damping provides a noticeable reduction 
of the structural vibration. 

In addition, when the electrohydraulic actuator is saturated, the movement of the 


end-point is slower than the ideal case without saturation. The simulation serves an 
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useful method to run the physical plant within the operating limit of the selected 
electrohydraulic actuator. 


B. RECOMMENDATIONS 

The long-term goal of this research is to design a controller being capable of 
achieving end effector positioning accuracy and stable movements. Simulation studies 
need to be continued for flexible multi- link system design using the ERLS model. In 
addition, the effects of the robustness and spillover problems need to be studied. 
Extention of the controller design and implementation to the multi-link of the ERLS 
model is eventually needed if the advantages of flexible manipulator is to be realized in 


practical industrial application. 
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APPENDIX A 


DERIVATION OF THE EQUATIONS OF MOTION FOR THE SINGLE - 
LINK FLEXIBLE MANIPULATOR 


The motion of the flexible link manipulator has been separated into large 
motions and small motions by the ERLS. The generalized coordinates are chosen to 
describe the large motions by the joint variable of ERLS and to describe the small 
motions bv nodal displacements of link. The two generalized coordinates chosen are 
the rigid body rotation,®, and the nodal displacement,U(0) which can be measured at 
the end of the link. The following are the two sets Lagrange equations used to develop 


the equations of motion; 
didt{ dKE/c0} - 0KE/o8 + OPE/0@ = F 


d/dt{ @KE/8U } - OKE/8U + OPE/GU = 0 


where 
KE - kinetic energy 
PE - potential energy 
§@- large motion joint variable measured between the ERLS link and 
the absolute x axis. 
LU - small motion displacement and slope 


F - generalized force for large motion 


The total kinetic energy is due to kinetic energy of the link motion (KE), and 
that of the actuators (KE),; 


KE =(KE), +(KE), 


The actuator are much more rigid than the link, the actuator is treated as a rigid 
body and all kinetic energy terms are calculated separatly. Generalized forces include 
any applied force and damping forces. The expressions utilized for the determination of 


the kinetic energy of the system are as follows; 


KE( arm) = 1/2 [ p R° {R} dv 


arm 
volume 
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KE( load) = 1/2 Trace J py Ry ¢ {Ry} dv 


load 
volume 


KE( rotor) = 1/2 Trace { He Re {R,} dv 


rotor 
volume 


The absolute position vector of the arm is determined from the following 
transformation; 


R= W(r+D) 


W = the 3 X 3 transformation matrix of function of theta,@, only. 
r 


the 3 X 1 local position vector of the arm measured from the 


coordinate system whose origin is the end of the ERLS link. 


d = the3 X | deformation vector that only includes the transverse 
displacement, derivation of this vector is discussed in chapter2. 
uu == the mass density of the steel arm. 


The absolute position vector of load is determined from the following 
transformation; 


Rr 
Dy 


the 3 X 3 transformation matrix due to the local deformation 
of the arm tip. 


ry = the 3 X 1 local position vector for the load 
ty = mass density of the steel load 


The absolute position vector of the hydraulic actuator is given by the following 
transformation; 


Rp =Apip 


Ap = the 3 X 3 transformation matrix due to the large motion rotation 
of the rotor 


rp = the 3 X | local position vector for the rotor 


Hp = the mass density of the actuator rotor, aluminum 
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Potential energy of the flexible arm comes from strain energy (PE), of the 
Structure due to deformation and the gravitational energy (PE), of the system. The 
Strain energy is produced by bending, tortion and extention. Along modeling, the 
lateral shear deflections are neglected so that the system energy of bending is due to 
flexure only. The expressions used for the determination of potential energy are as 
followes; 

PE, = i/2f EL, Wy ax 


link 
length 


t 
PES -fpr g dv 


link 
volume 


PEoL = - fy rh g dv 


load 
volume 
where 
E IL, = the flexural ngidity in the z direction, perpendicular to 
the plane of motion 
v= second derivative of the transverse displacement V with respect 


to the local x coordinate direction expressed interms of V(0) 


g = gravitational acceleration vector 


The following definitions for the inertia terms are utilized to simplify the 
computations and the resultant expressions in the equations motion; 


re = the3 X 3 inertia matnx of the load 
= See 
load 
volume 
Ip = the3 X 3 inertia matrix of the rotor 
= i) Up [Rr rR dv 
rotor 
volume 
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[L11(We' , Wg) = f wr’ Wo! Wo r dy 
link 
volume 
1112Wg',W) = f wrt We! WP dv 
link 
volume 
H21(W', Wg) = fp 8 We! Ww dv 
link 
volume 
H22.W',W) = fp Po wi w dv 
link 
volume 
122.Wt, Wp) = fa ¥o Wo We ¥ dv 
link 
volume 
1122(Wg', W) = fn Bo Wo! WY dv 
link 
volume 
1122(We', We) = fw '¥' Wo! Wo ¥ dv 
link 
volume 
ze 2 
= 2 
Pe ay 
Munere 
eine sem link Shidpe matrix 
Weg = derivative of the 3 x 3 matrix W with respect to the joint variable @ 


Wp = result of the simplification of the second time derivative of 


the transformation matrix W and is termed the residual acceleration. 


The following definitions are utilized to simplify the computations and the 


expressions used for computer coding of the equations of motion; 
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link 
length 


Hy = fprtdy 


link 
volume 


Ho, = Ju PE dv 


link 
volume 


as t 
Ha) = J Bp rps av 
link 
volume 
[’ = the second derivative of the shape function matrix 


C = the3 x 3 flexural rigidity matrix including only E 1¢2 


The actual derivation of the equations of motion requires detailed, tedious 
substitution of the expressions for the kinetic energy and potential energy into the 
Lagrange equations, and the Lagrangian formulation yields two sets of equations, One 
set describes the large motions and the other set describes the small motions. These 
two sets of equations are non- linear, coupled, second - order, ordinary differential 


equations represented as follows; 


Mog 9 + Mont SE OR OUE 


M qs echinacea alae 


where 

Mag = 1 X 1 effective inertia matrix for large motions 

M qn = 1 X 2 coupled inertia matrix of the small motion effect on large 
motions 

Ming = 2 X | coupled inertia matrix of the large motion effect on small 
motions 

Man = 2 * 2 effective inertia matrix for small motions 

Kn = 2 X 2 stiffness matrix for small motions 

Fg = 1 X | load vector for large motions 

F, = 2 X 1 load vector for small motions 
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8 = generalized coordinate of the large motions 
U = 2 X 1 generalized coordinate of the deformations representing 


the small motions-the displacement and slope. 
The coefficients for the terms are defined as follows: 


= t rC t r 
Mag = HII(Wg', Wg) + Ul 1122(Wg' , Wa) U 


i linaee ( Wo Dy Ip D,' We! ite AR Ip Ap’) 


= t 
Mon = Me Gvan a) ment, Lt M.),(M, L + Lyx thy) 


ie celeste rw + Hy) We + 
U'H,, Wg'g- Trace( Wg Dy; I, Dr‘ Ww, 
+ 2 Wg D, 1, Dy ' WwW 
ee t {t 


Mng = ( Trace ( Wg Dy Ty Dy yy" W"), 
Trace( Wg Dy Ip Dy yo! WS + 1121 WS,+1121(W", Wo ) 


ieee enemy) HiMya( v0) ) + (1,, + ly) * (9)) 


el t 
K, = Kitt + 1122( Wt, Wp ) 


= t t t yt 
ea H+, W LY ate W Dy, f, Dp yy W ), 
t tyt 
eee ecm i2) YS 


where 
L = the arm length 
My = the mass of the load 
M, = the first moment of the load with respect to the local y axis. 
T = the externally applied torque 


ArR = aresult of simplification of the second time derivative of the 
transformation matrix Ap and is termed D> residual acceleration. 
Ap® = a derivative of Ap with respect to the joint variable theta. 
Dy 1)PL12 = the derivatives of the arm tip deformation transformation matrix 
differentiated with respect to nodal deflection displacement and the nodal slope 


displacement respectively. 
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These expressions are coded and solved for in the computer program listed in 


Appendix B. 
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HI TLE 
CONST 


APPENDIX B 
SIMULATION PROGRAM 


ie Seo eZeel 


Kg 
GKV=2.0, GKP= 


) 
1.0 ,DEA=1.0 


Eeihaols COPY 


2 OR RR RR Rb ERR Rb RE ROR ROR Eb ER ob ER Ob bb FR RR Ob 2b OR Ob ob oF OF FR ob OF OF OF OF oF OF oF OF 


SIMULATION FOR DESIGN A CONTROL ALGORITHM 


THIS PROGRAM IS TO CALCULATE THE REQUIRED TORQUE AND REQUIRED 
CURRENT TO CONTROL THE END POSITION OF THE SINGLE- LINK ARM. 

THE ARM PARAMETERS ARE INPUTTED AND THE HYDRAULIC ACTUATION 
DYNAMICS ARE INCLUDED IN THE SIMULATION. 

[pee tNPUTS OR TRE SS rst&M ARE THE DESIRED ANGLE OF THE END OF ARM. 
THE OUTPUTS OF THE SYSTEM ARE THE TOTAL ANGLE, TOTAL ANGULAR 
VELOCITY, ANGULAR ACCELERATION WHICH ARE FED BACK TO THE 
COUTRCELER. THE CODINGSGONSISTS OF A MAIN 

PROGRAM AND SUBROUTINES AND ARE DESCRIBED BELOW. 


THE FOLLOWING PARAMETERS ARE DEFINED: 

Peach PPECTIVE CROSS-SECTIONAL AREA OF FLEXIBLE ARM 

Z.ARRDD<-3X3 SECOND TIME DERIVATIVE OF ROTOR RESIDUAL ACCELERATION 
MATRIX 

eee ers ROTOR TRANSFORMATION MATRIX DIFFERENTIATED WITH RESPECT 

Soe er reGiivE BULK MODULUS OF FLUID 

9. BIGF-2X1 RIGHT-HAND SIDE VECTOR FOR LARGE AND SMALL MOTION 
ACCELERATIONS 

6.BIGM-2X2 MATRIX OF LARGE AND SMALL MOTION ACCELERATION 
COR PRECIENTS 

(oe Pi-TOTALSEBACAGE COEFEICIENT OF THE ACTUATOR 

8.DEFM-DISPLACEMENT DEFORMATION VARIABLE 

> DEBMDSUINE DERIVATIVE OF DISPLACEMENT DEFORMATION VARIABLE 

10.DIFF,QOERR,QERR1,FACTOR-DUMMY VARIABLES 

11.DL1-343 DEFORMATION MATRIX 

P20 bel = 3s noe ORM nOheiea Chae DUPFERENTIATED WITH RESPECT TO THE 
DISPLACEMENT DEFORMATION VARIABLE 

13. DLID-=3xX3 FIRST TIME DERIVATIVE OF DEFORMATION MATRIX 

14.DM-ACTUATOR DISPLACEMENT 

Loe BoNODULUSPOF seunottecriy OF STEEL 

16.FN- RIGHT-HAND SIDE FOR SMALL MOTION ACCELERATIONS 

17.FQ-RIGHT-HAND SIDE FOR LARGE MOTION ACCELERATIONS 

18.G-3X1 GRAVITATIONAL ACCELERATION VECTOR 

fees A Se LiKe LenOue Nor INERTIA VECIOR 

Z0e nels LINKeoAreeMaAtR is FIRST MOMENT OF INERTIAVECTOR 

Zieooi~ nS LOADS RST MOMENT CF INERTIA VECTOR 

22e ners loOlau FELON PRESoURE COBFFICIENT 

Z23.PL-LOAD HYDRAULIC PRESSURE DROP 

2¢ Fo My DRAULTG SUPPLY PRESSURE 

Z25.QL-FLOW DELIVERED FROM THE SERVOVALVE 

26.SOL-2X1 VECTOR OF LARGE AND SMALL MOTION ACCELERATIONS 

Zi te tOROUE EPP PemeNGy 

28.TH-LARGE MOTION POSITION VARIABLE 

29.THD-TIME DERIVATIVE OF LARGE MOTION VARIABLE 

30.TORQUE-APPLIED TORQUE BY ACTUATOR 

31.U-241 ARM TIP DEFORMATION OF DISPLACEMENT 

32.UD-ARM TIP DEFORMATION DIFFERENTIATED WITH RESPECT TO TIME 

33.VT-TOTAL COMPRESSED VOLUME INCLUDING ACTUATOR LINES AND CHAMBERS 

34.W-3X3 LINK TRANSFORMATION MATRIX 

35.WD-3X3 FIRST TIME DERIVATIVE OF LINK TRANSFORMATION MATRIX 

36.WRDD-3X3 SECOND TIME DERIVATIVE OF LINK RESIDUAL ACCELERATION 
MATRIX 
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> Se a ea 


VVUVUVVUOVVVUOO 4 4H HHA HHH HAHA OH OOO 


37 ae TRANSFORMATION MATRIX DIFFERENTIATED WITH RESPECT TO 
oO ae eee FRACTIONAL AMOUNT OF INPUT CURRENT TO SERVO- 


39.XIINP-CURRENT INPUT EQUAL TO INITIAL AND FRACTIONAL AMOUNTS 
40.XIL-3X3 INERTIA MATRIX OF THE LOAD 
41.XIO-INITIAL INPUT CURRENT TO SERVOVALVE 
42.XIR-3X3 ROTOR INERTIA MATRIX 
43.XK11- PARTIAL LINK STIFFNESS 
44.XKN- LINK STIFFNESS 
45.XKV-SERVOVALVE SIZING CONSTANT 
46.XLL-LENGTH OF FLEXIBLE ARM 
47.XML-MASS OF LOAD 
48.XMNN-COEFFICIENT OF SMALL MOTION ACCELERATIONS IN THE 
SMALL MOTION DYNAMIC EQUATIONS 
49.XMNQ- COEFFICIENT OF LARGE MOTION ACCELERATIONS IN THE 
SMALL MOTION DYNAMIC EQUATIONS 
90.XMON- COEFFICIENT OF SMALL MOTION ACCELERATIONS IN THE 
LARGE MOTION DYNAMICS EQUATION 
a aoe COEFFICIENT OF LARGE MOTION ACCELERATION IN THE LARGE MOTION 
DYNAMICS EQUATION 
De Ee DUMMY MATRIX FOR USE IN FORMULATING THE EQUATIONS OF 
3.XMR-MASS OF ACTUATOR ROTOR 
54.XMU-MASS DENSITY OF STEEL FLEXIBLE ARM 
95.XMX-FIRST MOMENT OF LOAD WITH RESPECT TO THE LOCAL COORDINATE 
Y AAagS 
56.ZI-AREA MOMENT OF INERTIA OF FLEXIBLE ARM 
57.TAg- TOTAL total angle 
58.TAG1- TOTAL ANGULAR VELOCITY 
599.TAG2Z- TOTAL ANGULAR ACCELERATION 
60.RETO- REQUIRED TORQUE 
61.REIC -~REQUIRED CURRENT 





NITIAL 
INITIAL VALUES OF PARAMETERS ARE INPUTTED VIA XINIT SUBROUTINE 
DIMENSION U( 1) 399 (1), x90 tae ,DL1(3,3),WTH(3,3),ARTH(3,3), 
#XIR(3,3),XMON(1 ),UD G(3, 1) H21 (oe 
HWRDD(3, Orb (3,95 ,WDC3, 35 CO ARRDDY 3) Hal (1 3) RK), 
#XMNQ(1 ),W(3,3),XMNN(1 ),XKN(1 +) J FN(1 BIGM(2,2), 
H#BIGF(2,1),XIL(3,3) ,DL11(3,3) ,DEFMD(1) ,SOL(2), THD iy 
#n(1),E(1),Z1(1), FQ(1),GPOS(3) ,XITH(1), 
RAMU (L) ,XLL(1) , KML (1) AMR (1) XMK(1)/ TH(1) TORQUE (1), DEFM(1), 
gPS(1),CTM(1) VT (1) BE(1),DM(1),XKV(1),TE(1), 
HOL(1 PL(L) DEER CL) EINE(1) QGRRE 1), QERR(1),FACTOR(1), 
H#rEIC(1), rETO(1),FTT(1),CTB(1),FCO(1), 
#DEQ(1),DI(1),DEPL(1) 
FIXED I 
NOSORT 
D COMMON/FCDATA/C1,C2,  AlP,A2P,BETA1,BETA2 
C1=0.515462194 
C2=-0.205906794 
A1P=1.362220557 
A2P=0.981867539 
BETA1=1 .877920950 
BETA2=4.701142847 
CALL XINIT(TH,THD,DEFM,DEFMD ,VO,POSO,A,XML,XMU,... 
XLL,XMR,E,Z2I,PS ,CTM,VI,BE,DM,XKV,TE,QL,... 
PL,PLIC, REIC) 
DERIVATIVE 
* COEFFICIENTS FOR BOTH LARGE AND SMALL MOTION ACCELERATIONS 
* AND THE RIGHT-HAND SIDES ARE COMPUTED IN THE FOLLOWING 
‘ SUBROUTINES. ALSO,THE HYDRAULIC DYNAMICS ARE INCLUDED 


IN THE MAIN PROGRAM. 


a2 


x 
NOSORT 
k 


A HYDRAULIC DYNAMICS 
ae ne 
ie PE ech rs) GO TO 2 


3 OERRI (1 
OERR(1 
DIFF(1 
FACTOR(1 =i 
DIFF1(1)=DIFF( 


SORT 
PLI=INTGRL(PLIC, DIFF1,1) 
NOSORT 
PL(1)=PL1(1)/FACTOR(1 
TORQUE (1)=TE(1)*PL(1)*DM(1) 
MATRIX AND VECTOR FORMULATION SUBROUTINE 
CALL FORM(W,WTH,WD,DL1,DL1D,XIL,XIR,ARTH,WRDD,ARRDD,U,UD,... 
KUO PeereddyH21,DL11 .H41>XK11,A,KMU,XML,XLL,TH,THD,... 
,DEFMD ,E,Z1I,XMR,XMX_ ) 


COEFFICIENT OF LARGE MOTION ACCELERATION IN LARGE MOTION DYNAMICS 
EQUATION SUBROUTINE 


CALL XLMMQQ(XMQQ,U,XMQQP,DL1,WTH,ARTH, XIL,XIR,A,XMU,TH,DEFM ) 


COEFFICIENTS OF SMALL MOTION ACCELERATIONS IN LARGE MOTION DYNAMICS 
EQUATION SUBROUTINE 


CALL XLMMQN(XMOQN,A,XMU, XML, XLL, XMX, DEFM) 

RIGHT-HAND SIDE FOR LARGE MOTION DYNAMICS EQUATION SUBROUTINE 
CALL XLMFQ(FQ,U Ng P,DL1,WTH,ARTH, XIL,XIR,UD,H11,G, H21, WRDD, 
DL1D,WD,ARRDD,H41,TH,THD,DEFM,DEFMD ,A,XMU,XML,XLL,. 

TORQUE , FTT) 

LINK STIFFNESS MATRIX SUBROUTINE 

CALL SMKN(XKN,XK11,XMQQP,A,XMU, THD) 


COEFFICIENTS OF LARGE MOTION ACCELERATION IN SMALL MOTION 
DYNAMICS EQUATIONS SUBROUTINE 


ELL) SMMNQ(XMNQ ,DL1,WTH,XIL,DL11 7A, Layee, A,XMU,... 
x 


Nee! Ee Te) 








XKV ( 
a 
PO D@ er (1) -CIM(1)) 


* 
~ 


A Ot a Ob OF a a AO Ob OF a OF A: 


+ ob OF 


RIGHT-HAND SIDE OF SMALL MOTION DYNAMICS EQUATIONS SUBROUTINE 


CALL SMFN(FN,H21,W,G,WRDD,DL1,XIL,DL11, WE, PLID Hal Te, ... 
THD , DEFM, DEFMD) 


COEFFICIENTS OF SMALL MOTION ACCELERATIONS IN SMALL MOTION DYNAMICS 
EQUATIONS SUBROUTINE 


CALL SMMNN(XMNN, XMQQP,XML,A, XMU) 


ACCELERATION COEFFICIENTS MATRIX AND RIGHT-HAND SIDE VECTOR 
FORMULATION SUBROUTINE 


CALL BIGFOR(BIGM,BIGF,XMQQ, XMON, FQ, XMNQ, XMNN, XKN,FN,U,DEFMD) 
LINEAR EQUATION SOLVER FOR ACCELERATIONS SUBROUTINE 
CALL XLEQ(BIGM,BIGF,SOL) 


+ A Ft A 


A > 
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INTEGRATE ACCELERATIONS AND THEN VELOCITY TO GET LARGE MOTION 
ANGULAR POSITION AND SMALL MOTION,LOCAL COORDINATE,TIP POSITION 


+ HF 


DO 5 I=1,2 
SOLi (1)= =OL(1) 
5 CONTINUE 
SORT 
VEL=INTGRL(VO,SOL1,2) 
NOSORT 


THD(1)=VEL(1) 

DEFMD(1)=VEL(2) 

DOr LO fa 

VEL1(1)=VEL(I) 
10 CONTINUE 


SORT 

POS=INTGRL(POSO,VEL1,2) 
NOSORT 

TH(1)=P0S(1) 

ANG=TH(1) 

CUR=XIINP(1) 
DEFM(1)=P0S(2) 
ele oe 
AC2=SOL1(2 
VE1L=VEL(1 
VE2=VEL (2 
PO1=P0S(1 
PO2=P0S(2 


ANGULAR ACCELERATION ,VELOCITY ,ANGLE. 
TAG2=AC1+AC2/0.9985D0 
TAG1=VE1+VE2/0.9985D0 
TAG =PO1l +P02/0.9985D0 


* COEFFICIENT OF THE CONTROL ALGORITHM 
CALL CTE EGGe o eteT eG sade en 


+ OO A 


i DEFM,CTB,FCO, DEFMD 

; CALCULATE THE REQUIRED TORQUE 

I CALL REQTO(RETO,CTB,TAG1,TAG,FCO,GKV,GKP , DEA) 

: INVERSE HYDRAUIC DYNAMIC 

‘ CALCULATE THE REQUIRED CURRENT 
CALL REC(RETO,PL,TE,DM,VT, THD, PS ,CTM, XKV, REIC, BE) 
TORK=DETO(1) 
CU=DEner) 

TERMINAL 


METHOD ADAMS 

PRINT 12,0E=3) [aGe, lOnk ecu 
SAVE 6.0E-3,TAG 

os FINTIM=5.000 ,DELT=0.006 


PARAM 
END 

PARAM 
GRAPH 
LABEL 
GRAPH 
LABEL 
GRAPH 
LABEL 
GRAPH 
TABEL 
GRAPH 
LABEL 


GKV=4.0 ,GKP=4.0 


GKV=8.0,GKP=16.0 

G1,DE= TEK618) TIME (UN=SEC) , TORK 
G1) LOAD=2.1151 

G2,DE= TERGLB) TIME(UN=SEC), TAG 
62) TOAD=25 10S leKe 

G3, DE=TEK618)TIME, CU 

G3) LOAD=2.1151 KG 

G4 DE=TEK618) TIME,TAG1 

G4) LOAD=2.1151 KG 

gs DE=TEK618) TIME, TAG2 

G5) LOAD=2.1151 KG 
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END 
STOP 


a ElelInNG OF SUBROUTINES 


x 


FORTRAN 


* 
x 
SUBROUTINE XINIT(TH,THD,DEFM,DEFMD ,VO,POSO,A,ML,MU,LL, 
Toe ee omen vale bE OM, KY,16,0L,PL,PLIC ,REIC ) 
REAL*3 VO(2),POSO(2) ,ML,MU,LL,MR,TH,THD,DEFM,DEFMD, 
[ho oi Shon ome in Be.,DM,KV,TE,QOL,PL,PLIC, 
Tel C2 Ar Leora Degree wBeraAc,rEIC ,; 
TORE L, DEO: Un VE beep ere ll DETCZ 


ITORO= 44.554950510000000 
DM=6.2271D-05 
TE=.90000000000000 
PL=ITORO/ (DMATE) 
PLIC=PL 
CTM=3 .7064772D-13 
L=CTM*PL 
V=2.402963D-09 
PS=1.37888D+07 
VT=3 .05127D-04 
BE=690.D6 
A=6.17795D-04 
ML=2.1151000000000 
MU=7861.05000000000000 
LL=0.99850000000000 
MR=9.00011451000000 
E=2.0D11 
ZI=4.065D-10 
ee 0990000000 
V0(2)=0.000000000000000 
ee ee 
POSO(2)=-.14606414900000 
TH=P0S0(1) 
THD=VO(1) 
DEFM=P0S0 ( 
DEFMD=Vo(2 
REIC=4.0 
RETURN 
END 


2) 
) 


DOUBLE PRECISION FUNCTION ONE(X) 
REAL*8 C1,C2, A1P,A2P,BETA1, BETA2 
COMMON/FCDATA/C1 ,C2 ,A1P,A2P ,BETA1 , BETA2 


*eee¢¢ ¢ ¢ 


ONE= Sl vAlPs( C@c( BETAI*x)+COSH 


@2~(AZP~( COS (BETAZ*xX)+COSH( BETA2*X 
SIN(BETA2*X )+SINH( BETA2*X 


tt te tt 


EGRA (AzeH {0S BETAS CH) GOSH {BETA x 


RETURN 
END 


DOUBLE PRECISION FUNCTION TWO(X) 
REAL*S Cl C2, A1P,A2P,BETA1,BETAZ 
COMMON/FCDATA/C1,C2, A1P,A2P,BETAI, BETAZ2 


)+ 
+ 


) 


TWO= C1*(A1P* (COS (BETAL*X)+COSH 


# C2*(A2P* (COS ( BETA2*X ) +COSH( BETA2*X 
SIN( BETA2*X)+SINH( BETA2*X 
RETURN 

END 

DOUBLE PRECISION FUNCTION SIX(X) 

REAea omer ,C2y AlP ,A2P,BETA1, BETAZ 
COMMON/FCDATA/C1,C2, A1lP,A2P,BETA1, BETAZ2 
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( 
OECRERDBH Cos {eta ze *COsH | BETR2 X 





)+ 
ae 
) )*x2 


SIX= 985*(C1*(A1P* (COS (BETA1*X ee 











+ (BETAL "5 J SIN(BETA1*X)+SINH(BETA1 + 
# COR (naps CoS (BETA2*x )+COSH BETAD*S 
+ SIN BETAZS +SINH(BETA2*X ye 
i Xx (C1*(A1P* (COS (BETAL*X +COSH 
+ (BETAL*X) )+ SIN(BETA1*X)+SINH(BETAL + 
# * (A2P* (COS (BETA2*X ) +COSH BETA2*R 
4 SIN(BETA2*K)+SINH(BETA2*X ) 

RETURN 

END 

DOUBLE PRECISION FUNCTION EIGHT(X) 

REAL*8 C1,C2, AlP,A2P,BETAI1 , BETA2 

COMMON/FCDATA/C1,C2, A1P,A2P,BETA1, BETA2 

EIGHT= C1*BETA1*BETA1*(A1P*(-COS 

(BETA1*X)+COSH(BETA1*X) )+ SIN) BEALS pS NH BETA Ef 
+ C2*BETA2*BETA2* (A2P*( -COS (BETA2*X) +COSH(BETA2*X 
(-SIN(BETA2*X)+SINH(BETA2*x) ) ) 
RETURN 
END 


SUBROUTINE FORM(W,WTH,WD,DL1,DL1D,XIL,XIR,ARTH,WRDD,ARRDD,U,UD, 


#HIQOP, ne sell .H21, DL1 I. /H41,XK11,A,MU,ML,LL, TH, THD ,DEFM,DEFMD, 
REAL*8 W(3, a RS 3) ab AD (3 , 2), Dias .3) sDLID(3, 3) ,XIL(3,3) 
REAL*8 hae /ARTH ( D(3,3),A 3,3).U UD 
REAL*8 XM HOG S3, 1) ee 35, HoT CL. 03) DL11(3,3) 

REAL*8 H41(1, i. XK11 /MU,ML,LL,MR,MX, TH 
REAL*8 THD, DEFM, DEFMD Aenea 


REAL*8 C1,Al1P,BETAI 
REAL*8 ONE,TWO,EIGHT,Y 
EXTERNAL ONE, TWO ,EIGHT 


W(1,1)=1.00000000000000 
W(1,2)=0.00000000000000 
W(1,3)=0.00000000000000 
W(2/1)=LL*DCOS (TH) 
W(2,2)=DCOS(TH) 
W(2,3)=-DSIN(TH) 
W(3,1)=LL*DSIN(TH) 
W(3,2)=DSIN(TH 
W(3,3)=DCOS(TH 
WTH(1,1)=0.00000000000000 
WTH(1,2)=0.00000000000000 
WTH(1,3)=0.00000000000000 
WTH(2,1)=-LL*DSIN(TH) 
WTH(2,2)=-DSIN TH 
WTH(2,3)=-DCOS 

WTH(3,1 =LLeDeOS (TH) 
WTH(3,2)=DCOS(TH) 
WTH(3,3)=-DSIN(TH) 
WD(1,1)=0.00000000000000 
WD(1,2)=0.00000000000000 
WD(1,3)=0.00000000000000 
WD(2.1)=-LL*DSIN(TH)*THD 
WD(2,2)=-DSIN(TH) *THD 
WD(2,3)=-DCOS (TH) *THD 
WD(3,1)=LL*DCOS (TH) *THD 
WD(3,2)=DCOS(TH)*THD 
WD(3,3)=-DSIN(TH) THD 
DL1(1,1)=1.00000000000000 
DLI@.2 = 00000000000000 
DL1(1,3)=0.00000000000000 
DL1(2,1)=0.00000000000000 
DL1(2,2)=1.00000000000000 
DL1(2,3)=0.0000000000000 
DL1(3,1)=DEFM 
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Se aT ee aerate 
DL1(3,3)=1.00000000000000 
DL1D(1,1)=0.00000000000000 
DL1D(1,2)=0.00000000000000 
DL1D(1,3)=0.00000000000000 
DL1D(2,1)=0.00000000000000 
DL1D(2,2)=0.0000000000000 
DL1D(2,3)=0.000000000000 
DL1D(3,1)=DEFMD 
DL1D(3,2)=0.000000000000 
DL1D{3,3)=0.09000000000000 
VI ieee L 

XIL(1,2)=0. 01678545900000 
XIL(1,3)=.00000000000000 

WE 2 ou =0 .01678545900000 
XIL(2,2)=1.77679D-04 

ME 23 =0 .00000000000000 
XIL(3,1)=.00000000000000 
Xie 2 =0.60000000000000 
XIL(3,3)=2.986791D-03 
MX=.01678545900000 
XIR(1,1)=MR 
XIR(1,2)=0.000000000000000 
XIR(1,3)=0.000000000000000 
XIR(2,1)=0.000000000000000 
XIR(2,2)=.027467130000000 
PIR( 253 =0 .000000000000000 
XIR(3,1)=0.000000000000000 
XIR(3,2)=0.000000000000000 
XIR(3,3)=.02746713000000 
ARTH(1,1)=0.00000000000000 
ARTH(1,2)=0.00000000000000 
ARTH(1,3)=0.00000000000000 
ARTH(2,1)=0.00000000000000 
ARTH(2,2)=-DSIN(TH 
ARTH(2,3)=-DCOS(TH 
ARTH(3,1)=0.00000000000000 
ARTH(3,2)=DCOS(TH) 
ARTH(3,3)=-DSIN(TH) 
WRDD(1,1)=0.00000000000000 
WRDD(1,2)=0.00000000000000 
WRDD(1,3)=0.00000000000000 
WRDD(2,1 =-LL*DCOS ( eee 2) 
WRDD(2,2 =-DCOS (TH) * ( THD*%2) 
WRDD(2,3)=DSIN(TH)*(THD**2) 
WRDD(3,1 =-LL*DSIN(TH)*(THD**2) 
Wee es ere cans eee 
WRDD(3,3)=-DCOS(TH)*(THD**2 
ARRDD(1,1)=0.00000000000000 
ARRDD(1,2)=0.00000000000000 
ARRDD(1,3)=0.00000000000000 
ARRDD(2,1)=0.00000000000000 
ARRDD(2,2 =“Dcos (TH) (THD*2) 
ARRDD(2,3)=DSIN(TH) *(THD**2) 
ARRDD(3,1)=0. 6099000000000¢ 
ARRDD(3,2 =<DSIN( TH)’ es 
ARRDD(3,3)=-DCOS(TH) *(THD**2 
U =DEFM 

UD =DEFMD 


CALL DQG4(-LL,0.DO,TWO, Y) 
XMOOP =Y 
ete {\=9;89099009000099 
D 


=4.85651900000000 
=-2.42825869300000 


=0.00000000000000 
ie 1 
it Za 
Hl 3 =0.00000000000000 
2 
Q 


-9.80660000000000 
2 


eo 
CALL 


3, 
i 
1 
1 
1 =0.00000000000000 


=0.00000000000000 
G4(-LL,-0.D0,ONE ,Y) 


1 
I. 
t 
i: 
J 
i 


ey 


+ + + 


+ A + 


H21(1,3)=A5M MU*Y 
0 I=1,3 
DO 60 J=1 3 
DL11(I,J)=0.00000000000000 
60 CONTINUE 
50 CONTINUE 
DL11(3, 1)= 1.00000000000000 
H41(1,1)=ML 
H41(1,2 oF 000000000000000 
H41(1.3)=.016785459000000000 
CALL DQG4(-LL,0.DO,EIGHT,Y) 
XX11 =Y*ERZI 
RETURN 
END 


SUBROUTINE ra L1,WTH,ARTH,X 
REAL*8 beta Ur pened 3,3 3) / ARTE 
REAL*8 P2(3,3),P3(3 P5(3,3),P6(3, 
REAL*8 U Hage roid 3), ee 3), 
REAL*8 XIR(3,3),A,T P 

M=2 

tel 

N=3 

Moos =0.00000000000000 

CALL TRANS(U,UT,L,L) 

CALL ea a /XMOOP,L,L,L,P) 

CALL MATMUL(P,U,L.L,L,SP) 

CALL TRANS(DL1,DL1T,N,N) 
CALL TRANS(ARTH,ARTHT,N,N) 
CALL TRANS(WTH ,WTHT,N,N) 
CALL MATMUL(WTH,DL1,N,N,N 
CALL MATMUL(P1,XIL,N 
CALL MATMUL(P2,DL1T, 
CALL MATMUL(P3,WTHT, 
CALL MATMUL(ARTH, AIR 


~ WWRNO 


H 





CALL MATMUL(P5,ARTH 
CALL MATADD P4,P6,N 


END 


SUBROUTINE ADEMON ANON 2 ,A,MU,ML,LL,MX ,DEFM ) 
REAL*8 XMON DEFM 

REAL*8 Cl,A1P,BETA1,Y, aoe ,C2,A2P ,BETA2 

EXTERNAL SIX 


CALL DOG4(-LL,0.D0,SIX ¥) 
XMON SCAMHO RE) 4 CMLALL +MX 
RETURN 

END 


SUBROUTINE XLMFO(FQ,U,XMQQP,DL1,WTH,ARTH,XIL,XIR,UD,H11,G,H21, 
#WRDD,DL1D,WD,ARRDD, H41, It THD, DEFM,DEFMD ,A, MU, ML , nie 
#TORQUE, /FTT) 
REAL*8 


Py P. ‘Pa 1,3),P2(1,3),P3(1,3),P4(3,3) ,P5(3,3) 
REAL*8 P P7(3,3),P8(3,3),P9(3,3),P10(3,3),P11(1,3),P12(1,3) 
REAL*8 saa Ae T(3,3),DLIDT(3,3),WDT(3,3) ,ARRDDT(3, 3) 
REAL*8 UT,DL1T(3,3),WRDDT(3,3),FPF(3,3),FPS(3,3) ,WTHT(3,3) 
REAL*8 U NQOR -DL1(3,3),WTH(3,3) ARTH(3, 3) XEL(G,3) 
REAL*8 zip 3) (Ub 2) H21(1,3).WRDD(3, 3) 
REAL*8 DL1D(3,35,WD(3,3) ,-ARRDD(3,3),H41(1,3),MU,LL,ML 
REAL*8 A, TORQUE,FQ,TH, THD, DEFM, DEFMD 
REAL*8 FP,SP,TP.TEP,FTHP ,FTT 


M=2 


a > 


L=1 
N=3 
CALL TRANS (UU ieee) 
FOP= KMOOP 
CALL MA HUL (UT, /FOP,L, 
CALL MATMUL(P,UD,L,L, 
CALL TRANS(WTH,WTHT,N, 
CALL MATMUL(H11,WTHT,L, 
CALL MATMUL(P1,G,L,N.L, 
CALL MATMUL(UT.H21,L.M_N, 
aL MATMUL (P2 WIHT L NON 
CALL MATMUL(P3,G,L,N, 
Gane TRANS (DLI DLT N 
CALL TRANS (WRDD,WRDDT 
CALL MATMUL(WTH,DL1,N 
CALL MATMUL(P4,XIL,N,N.N, 
CALL MATMUL(PS,DLIT,N,N,N 
CALL MATMUL{P6,WRDDT,N 
CALL TRANS (DL1D DLIDE 
CALL TRANS(WD,WDT,N,N) 
CALL MATMUL(WTH,DL1,N, 
CALL MATMUL(P7,XIL,N,N 
N 
N 


Zoos tt 
ZZ Zz 





CALL MATMUL(P8,DL1DT, 
CALL MATMUL(P9,WDT,N, 
CALL TRANS(ARRDD,ARRDD 
CALL a ert xTR 
CALL MATMUL(P10,ARRDDT 
DO 30 I=1,3 
Domagns=1,3 
FPS(I,J)=FPS(I,J)*2. 


40 CONTINUE 
30 CONTINUE 


CALL HazapD (FE FPS,N,N, see 
CALL MATADD(FPFH, FPT,N,N, FHP 

CALL TRACE(FHP,N,TFP 

CALL MATMUL(H41, e231) 

CALL MATMUL(P11,WTHT,L,N,N,P12 

CALL MATMUL(P12,G,L,N,L,FTHP) 
FTT=(-2.*A*MU*FP) + SP + TP - TEP + FTHP 
FO=FTT+TORQUE 

RETURN 


END 





SUBROUTINE SMKN(XKN,XK11 MOOR, A,MU,THD) 
REAL*8 XKN,KNP atop mat, A, THD ,MU 
KNP =XMQQP ey AMU* (THD* 

XKN =KNP +XK 

RETURN 

END 


oa SMMNQ(XMNQ,DL1,WTH, XIL,DL11,W,TH,DEFM,A,MU, 


REAL*8 xa ICL Ie 2) AT (3, 3), /P1(3,3),P2(3,3) 
REAL*8 P3( 524 th 3 ae 3,3), WTH(3,3),XIL(3,3) 
REAL*8 W(3,3 ‘TH, DEFM ,A,MU,LL 
REAL*8 Cl, ae Se 

REAL*8 SIX,Y ,C2,A2P,BETA2 

EXTERNAL SIX 

M=2 

fea 

N=3 

Onis ae /DL11T,N,N) 

CALL TRANS(W,WT,N,N 

CALL MATMUL(WTH,DL1,N,N,N,P1) 

CALL MATMUL(P1,XIL,N,N,N,P2) 

CALL MATMUL(P2,DL11T,N,N,N,P3) 


719 





+ + 4 


re 


20 
10 


50 
10 


eine MATHUL (F3 WT,N,N,N,P4) 
ene TRACE (P4,N /TEP1)~ ; 
CALE DQG4 (- -LL,0.D0,SIX ,Y) 


ae FRleY (PRAM ) 
RETURN 
END 
SUBROUTINE SMEN(FN,H21,W/Ga@WRDD,DLI XIE, DLIIy WD ,DL1D,H41,TH, 


#THD , DEFM, DEFMD ) 


REAL*8 =N, Pies yeas, uy Bose 3) ,P4(3,3) 2s seeps (oS) 
REAL*8 P7(3,3) P38 (3, 3), Po (337 


REAL*8 Pigit eo) r ue /TP,FP ,SP 
REAL*8 FN1(3.3), i Haid, 3), WEDD(3, 3) DL IDK 363) 
REAL*8 WD(3,3),H41(1,3) ,XIL(3,3) ,W(3,3) ,DL11(3,3) 
REAL*8 DL1(3,3),DL117T(3,3),WT(3,3 
REAL*8 TH, THD ,DEFM, DEFMD, TEN1, FN3 
L=1 

N=3 

CALL TRANS (W,WT,N,N) 

earn MATHUL(H21 'WL,L,N,N,P1) 

CALL MATMUL(P1,G,L, NAL EP) 

CALL TRANS(DL11,DL11T,N,N 

CALL MATMUL(WRDD ,DL1,N,N,N,P2) 

CALL MATMUL(P2,XIL,N,N,N,P3) 

CALL MATMUL(P3,DL11T,N,N,N,P4) 

CALL MATMUL(P4,.WT,N,N,N,PS5) 

CALL MATMUL(WD ,DLID,N.N.N,P6) 

CALL MATMUL(P6,XIL,N,N,N,P7) 

CALL MATMUL(P7,DL11T,N,N,N,P8) 

CALL MATMUL(P8,WT,N,N,N,P9) 

DO 10 I= 1,3 

DO 20 J=1,3 

P9(I,J)= P9(I,J)*2. 

CONTINUE 

CONTINUE 

CALL MATADD(PS,P9,N,N,FN1L) 

CALL TRACE(FN1,N "EEN § 

SP =TFN1 

CALL MATMUL(H41,DL11T,L,N,N,P14) 
CALL MATMUL(P14,WT,L,N,N,P15) 

CALL MATMUL(P15,G,L,N,L,EN3) 

TP =FN3 

FN=FP -SP + TP 

RETURN 

END 

SUBROUTINE SMMNN(XMNN, XMQQP ,ML,A,MU ) 
REAL*8 XMNN,XMOOP ,ML,MU,A 

XMNN =ML 

XMNN =XMNN + XMQQP *A*MU 

RETURN 

END 

SUBROUTINE MATMUL(A,B,M,L,N,C) 
REAL*8 A(M,L),B(L,N),C(M,N) 

DO 10 I=1,M 

DO 20 J=1,N 

C(I,J)=0.0 

DO 30 INDEX=1,L 

C(I, J)=C(I,J) + A(I, INDEX) *B(INDEX, J) 
CONTINUE 

CONTINUE 
CONTINUE 
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A A 


a+ A + 


A A + 


a Ot 


a A 


Pe 


eoeeeeeoe*eee54s5usti@e#eeo#tee%##enee#e#eekhmeie#e#eeekeeeé 


On PON, INO , XMNN , XKN,FN, UEREEMD) 


é ? 


18) 


END 
SUBROUTINE TRANS(A,B,M,L) 
REAL*8 A(M,L),B(L,M) 
DO 10 I=1,M 
DO 20 J=1‘L 
B(J,I)=A(I,d) 

20 CONTINUE 

10  CCNTINUE 
RETURN 
END 
SUBROUTINE TRACE(A,M, TRAC) 
REAL*8 A(M,M) 
TRAC=0.0 
pO 10 I=1,M 
TRAC=TRAC + A(I,I) 

10 CONTINUE 
RETURN 
END 
MATRIX ADDITION SUBROUTINE 

SUBROUTINE MATADD(A,B,M,L,C) 
REAL*8 A(M,L),B(M,L),C(M,L) 
DO 10 I=1,M 
DO 20 J=1,L 
C(I,J)=A(I,J) + B(I,d) 

20 CONTINUE 

10 CONTINUE 
RETURN 
END 
SUBROUTINE BIGFOR(BIGM, arg Mo 
REAL*8 BIGM(2,2),BIGF(2,1 
REAL*8 FN MOQ, FO, bEEMD 
M=2 
L=1 
BIGM(1,1)=MQQ 
BIGM(1,2 =KNON 
BIGM(2,1)= 9 
BIGM(2,2)= 
BIGE (isl = FO 
CALL MATMUL(XKN,U,L,L,L,P) 
BIGF(2,1)=EFN =P 
RETURN 
END 
SUBROUTINE XLEQ(BIGM,BIGF,SOL) 
REAL*8 BIGM(2,2),BIGF(2,1),SOL(2) , WKAREA( 
N=2 
CALL LEQT2F (BIGM,M,N,N,BIGF,M,WKAREA, IER) 
DO 10 I=1,2 
SOL (1 )= =BIGE(T, 1) 

10 CONTIN 

RETURN 
END 


SUBROUTINE LEQT2F (IMSL) 


8] 


HHH KSA SEE ee aK OD OO OO OM RO OD OD OO OO OO OO OD A OO FO OE OC 


+ 


PURPOSE 


USAGE 
ARGUMENTS A 


Has 


EveT 


WKAREA 
DER 


REQD. IMSL ROUTINES 


LINEAR EQUATION SOLUTION - FULL STORAGE 
MODE - HIGH ACCURACY SOLUTION 


CALL LEQT2F (A,M,N,IA,B,IDGT,WKAREA, IER) 


INPUT MATRIX OF DIMENSION N BY N CONTAINING 
THE COEFFICIENT MATRIX OF THE EQUATION 


A : 
- NUMBER OF RIGHT-HAND SIDES. (INPUT) 


ORDER OF A AND NUMBER OF ROWS IN B. (INPUT) 


- ROW DIMENSION OF A AND B EXACTLY AS SPECIFIED 


IN THE DIMENSION STATEMENT IN THE CALLING 
PROGRAM. (INPUT 
INPUT MATRIX OF DIMENSION N BY M CONTAINING 
THE RIGHT-HAND SIDES OF THE EQUATION AX = B. 
ON OUTPUT, THE N BY M MATRIX OF SCLUTTGIE 
REPLACES B. 
INPUT OF TION: 
IF IDGI IS GREATER THAN 0, THE ELEMENT Seen 
A AND B ARE ASSUMED TO BE CORRECT TO IDGT 
DECIMAL DIGITS AND THE ROUTINE PERFORMS 
AN ACCURACY TEST. 
IF IDGT EQUALS 0, THE ACCURAG SEs? Is 
BYPAsSEL. 
ON OUTPUT, IDGT CONTAINS THE APPROXIMATE 
NUMBER OF DIGITS IN THE ANSWER WHICH 
WERE UNCHANGED AFTER IMPROVEMENT. 
WORK AREA OF DIMENSION GREATER THAN OR EQUAL 
TO N**2+3N. 
ERROR PARAMETER. (OUTPUT) 
WARNING ERROR 
IER = 34 INDICATES THAT THE ACCURACY TEST 
FAILED. THE COMPUTED SOLUTION MAY BE IN 
ERROR BY MORE THAN CAN BE ACCOUNTED FOR 
BY THE UNCERTAINTY OF THESDete. hte 
WARNING CAN BE PRODUCED ONLY IF IDGT [5 
GREATER THAN 0 ON INPUT. (SEE THE 
CHAPTER L PRELUDE FOR FURTHER DISCUSSION. ) 
TERMINAL ERROR 
IER = 129 INDICATES THAT THE MATRIX IS 
ALGORITHMICALLY SINGULAR. (SEE THE 
CHAPTER L PRELUDE). 
IER = 131 INDICATES THAT THE MATRIX IS TOO 
ILL-CONDITIONED FOR ITERATIVE IMPROVEMENT 
10 BE JERFECH Vee 
SINGLE/LUDATN , LUELMN , LUREFN, UERTST , UGETIO 
DOUBLE/LUDATIN , LUELMN, LUREFN, UERTST, UGETIO, 
VXADD , VAMUL , VXSTO 


SUBROUTINE LEQT2F (A,M,N,IA,B,IDGT,WKAREA, IER) 


DIMENSION 
DOUBLE PRECISION 


DO 5 ‘T=1,N 


A(IA,1),B(IA,1),WKAREA(1) 
AB WKAREA,D1, D2,WA 
FIRST EXECUTABLE STATEMENT 
INITIALIZE IER 


WKAREA (JJ) =A(I,L) 


JJ=HS I+ 
5 CONTINUE 


DECOMPOSE A 


82 


+ + 


+ oe Ee a a A ee es ee 


+ 


CALL LUDATN Need fi ge EAA I), WRAREA(K) , 


IF (TER. (e1e 128) GO TO 25 
IF (IDGT .EQ. 0 .OR. IER .NE. 0) KK = 1 


Doms f= 1° 
PERFORMS THE ELIMINATION PART OF 


AX = B 
CALE LUBEMN (A, IA,N,B(1,L), WKAREA(J), WKAREA(MM) ) 
REFINEMENT OF SOLUTION TO AX = B 


IF (KK .NE. 0) 
* CALL LUREFN (WKAREA,N,N,A,IA,B(1,1),IDGT,WKAREA(J) , WKAREA(MM), 
x oe WKAREA(K) , WKAREA(K) , JER) 

D =] 


WN 
B(II,I) = WKAREA(MM1+II) 
10 CONTINUE 
IF (JER.NE.0) GO TO 20 
15 CONTINUE 
Core 25 
20 IER = 131 


DO 30 I = 1, 
te J) =ifKBREA( JI) 
ie 

30 CONTINUE 
IF (IER .EQ. 0) GO TO 9005 
9000 CONTINUE 
CALL UERTST (IER,6HLEQT2F) 
9005 RETURN 
EN 


SUBROUTINE DQG4 (IMSL SUBROUTINE) 


PURPOSE 
TO COMPUTE INTEGRAL(FCT(X), SUMMED OVER X FROM XL TO XU) 


USAGE 

CALL DQG4 (XL,XU,FCT,Y) 

PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT 
DESCRIPTION OF PA ETERS 


XL - DOUBLE PRECISION LOWER BOUND OF THE INTERVAL. 
XU = S0OUBLE FRECISION UPPER BOUND OF THE INTERVAL. 
EGt - THE NAME OF AN EXTERNAL DOUBLE PRECISION FUNCTION 
SUBPROGRAM USED. 
iy - THE RESULTING DOUBLE PRECISION INTEGRAL VALUE. 
METHOD 


EVALUATION IS DONE BY MEANS OF 4-POINT GAUSS QUADRATURE 
FORMULA, WHICH INTEGRATES POLYNOMIALS UP TO DEGREE 7 
EXACTIEY SenOheReeerRenCE, SEE 

V.I.KRYLOV, APPROXIMATE CALCULATION OF INTEGRALS, 
MACMILLAN, NEW YORK/LONDON, 1962, PP.100-111 AND 337-340. 


SUBROUTINE DQG4(XL,XU,FCT,Y) 


Reeeas PRECIOLEON@S Gn: &,8,C,FCT 

- 5DO* (XU+XL) 
Ba XU- XL 
C=.43056815579702629D0*B 
Y=.17392742256872693D0*(FCT(A+C)+FCT(A-C) ) 
C=.16999052179242813D0*B 
Y=B* (Y+. 32607 257743127307D0* (FCT (AtC)+FCT(A-C))) 
as 


aDEFMD) CCO(XMOQ, KMON” XMNN] XMNOQ| FIT) FN) XKN,DEFM,CTB, FCO, 
M 
REAL*8 XMQQ,XMNN, XMQN, FIT, FN,XKN,U,CTB,FCO,CTA,B11,B22 


83 


READ“ 6gbo> eae oe F22,DEEMND Poo e 


B11=XM0Q/0 
B22=KM Ne ‘Bil 


F22= aa inae 


F33=(CTA*0.00)*DEFMD 
FCO=F1l1-F22-FTT -F33 
RETURN 

END 


SUBROUTINE See ene ,CTB, TAG1, TAG, FCO, GKV, GKP, DEA) 
REAL coms DETOZ2 rETO, GKV, GKP, DEA, THGIe TAG, FCO, CTs 
DETO1=(- GKVATAG1 )+GKP*(DEA- Tie) 

DETO2=CTB*DETO1L 
RETO=DETO2Z+FCO 
RETURN 

END 


SUBROUTINE REC(RETO,PL,TE,DM,VT,THD,PS,CTM,XKV,DEIC,BE) 
REAL*8 RETO,PL, TE,DM,VT,THD,PS,CTM,XKV,REIC ,BE 
REAL*8 DI,DEQ, DEFF DEICl, DEIC2 
DEPL =DETO te 
Di vs DEPL PL )/(4.0D0*BE  *0.006D0) 
DEQ =DM *THD +CTM *DEPL +DI 
P= (PS - DEPL ) 
IF(DEFP .LT.0eG))eGommo U3 
DEIC1=DSORT (DEFP) 
GO TO a4 
13 DEFP=- 
DEICL= =DSQRT(DEFP) 
14 DEIC2=DEIC 
REIC =DEQ /(DEIC2 ) 
RETURN 
END 
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